# -*- coding:utf-8 -*-
"""
Created on 2009-10-28

@author: Qiu wenfeng
"""
from neumannboundcond import *
from dirac import *
from curvature_central import *

import scipy
import scipy.signal


def lsecml05(u0, g, lmd, mu, alf, epsilon, delt, numIter):
    '''
    this function refers to ChunmingLi matlab code. I rewrite to python
    lsecml05(u0, g, lmd, mu, alf, epsilon, delt, numIter) updates the level set function
    according to the level set evolution equation in Chunming Li et al's paper:
        C. Li, C. Xu, C. Gui, and M. D. Fox "Level Set Evolution Without Re-initialization: A New Variational Formulation"
        in Proc. IEEE Conf. Computer Vision and Pattern Recognition (CVPR'05),
        vol. 1, pp. 430?36, 2005.
    Usage
    @param u0: level set function to be updated
    @param g: edget indicator function
    @param lambda: coefficient of the weighted length term L(\phi)
    @param mu: coefficient of the internal (penalizing) energy term P(\phi)
    @param alf: coefficient of the weighted area term A(\phi), choose smaller alf 
    @param epsilon: the papramater in the definition of smooth Dirac function, default value 1.5
    @param delt: time step of iteration, see the paper for the selection of time step and mu 
    @param numIter: number of iterations.    
    @copyright: Chunming Li, Qiu wenfeng
    @author: Qiu wenfeng    
    '''
    u=u0
    vy, vx = scipy.gradient(g)
    for k in range(numIter):
        u = neumannboundcond(u)
        uy,ux=scipy.gradient(u)
        normDu = scipy.sqrt(ux**2+uy**2+1e-10)
        Nx = ux/normDu
        Ny = uy/normDu
        diracU = dirac(u,epsilon)
        K=curvature_central(Nx,Ny)
        
        weightedLengthTerm=lmd*diracU*(vx*Nx + vy*Ny + g*K)
        
        laplacian = scipy.array([[0,1,0],[1,-4,1],[0,1,0]],float)
        deriv2 = scipy.signal.convolve2d(u,laplacian,mode='same',boundary='symm')
        penalizingTerm=mu*(deriv2-K)
#        penalizingTerm=mu*(4*del2(u)-K) # del2???

        weightedAreaTerm=alf*diracU*g
        # update the level set function
        u=u+delt*(weightedLengthTerm + weightedAreaTerm + penalizingTerm)
    return u